function [edge_new] = genimganddetectsubpixeledge(imageSize, imageEdge, fitMethod, nNeighborPts)
%GENIMGANDDETECTSUBPIXELEDGE 亚像素边缘提取算法，一维边缘提取仿真
%
% Input Arguments
% # imageSize：2*1矩阵，图像的高度m和宽度n
% # imageEdge：标量，图像边缘，图像的宽度为m中的一点为图片分界点
% # fitMethod：int8类型标量，表示曲线拟合算法，1为二阶多项式拟合，2为取对数拟合，3为高斯拟合
%
% Output Argumnets
% # edge_new：长度为m的列向量，图片亚像素边缘，理论值为Edge+0.5
% 
% Reference
% # 张虎, 达飞鹏, 邢德奎: 光学测量中椭圆圆心定位算法研究

imageOriginal = zeros(imageSize(1), imageSize(2)) + 50;
imageOriginal(:, imageEdge+1:imageSize(2)) = 150;
sigma = 2;
gaussFilter = fspecial('gaussian', [49, 49], sigma);
image_new = imfilter(imageOriginal, gaussFilter, 'replicate');
edges = edge(image_new, 'Canny');
B = bwboundaries(edges, 'noholes');
edge_new = NaN(size(B{1, 1}));

for iEdgePt=1:1 % size(B{1,1}, 1)
    thisEdgePt = double(B{1, 1}(iEdgePt, :)');
    ind = NaN(nNeighborPts*2+1, 1);
    grayScale = NaN(nNeighborPts*2+1, 1);
    ind(nNeighborPts+1, 1) = 0;
    grayScale(nNeighborPts+1, 1) = image_new(thisEdgePt(1, 1), thisEdgePt(2, 1));
    for iNb=[-nNeighborPts:-1 1:nNeighborPts]
        pNb = thisEdgePt + [0; iNb];
        ind(iNb+nNeighborPts+1, 1) = iNb;
        grayScale(iNb+nNeighborPts+1, 1) = image_new(pNb(1, 1), pNb(2, 1));
    end
    GSDiff = diff(grayScale);
    ind = ind(1:(end-1), 1) + 0.5;
    signIsDiffWithMajorSign = (sign(GSDiff)~=sign(grayScale(end, 1)-grayScale(1, 1)));
    GSDiff(signIsDiffWithMajorSign) = NaN;
    ind(signIsDiffWithMajorSign) = NaN;
    ind(isnan(GSDiff)) = []; % 去掉无效点
    GSDiff(isnan(GSDiff)) = [];
    
    if fitMethod == int8(1)
        c = intcoeffs(ind+0.5) - intcoeffs(ind-0.5);%二阶曲线拟合
        ABC = c \ (GSDiff);
        delta = -ABC(2, 1) / (2*ABC(1, 1)); % mu=-B/(2*A), sigma=sqrt(-1/(2*A))
    elseif fitMethod == int8(2)
        c = logarithm(ind);
        ABC = c \ log(GSDiff);
        delta = -ABC(2, 1) / (2*ABC(1, 1)); % mu=-B/(2*A), sigma=sqrt(-1/(2*A))
    elseif fitMethod == int8(3)
        %高斯曲线拟合
        p = lsqcurvefit(@gauss, [0 1 1], ind, GSDiff); % TODO: 预估p初始值
        delta = p(1) - 0.5;
    end
    edge_new(iEdgePt, :) = thisEdgePt' + [0, delta];   
end
end

function [coef] = intcoeffs(d)
coef = [d.^3/3, d.^2/2, d]; % A*d^2 + B*d + C的原函数为(A*d^3)/3 + (B*d^2)/2 + C*d
end

function [coef] = logarithm(d)
coef = [d.^2, d, d.^0];%采用取灰度值,不取积分，直接取delta处为灰度差值，即取中点值为单位长度的积分
end

function [y] = gauss(p, k)
ur = p(1);% 高斯曲线拟合
sigma = p(2);
a = p(3);
x = (k + 0.5 - ur) ./ sigma;
y = a * 1/sqrt(2*pi) * exp(-x.^2/2);
end